Introduction to geopandas and cartopy#

Basic Setup#

Again we will be using pandas and matplotlib.

import pandas as pd
import matplotlib.pyplot as plt

We’ll also supress a few disturbing warnings.

import warnings
warnings.filterwarnings('ignore')

Why do we need something other than pandas?#

Let’s reload again our example dataset of conventional power plants in Europe as a pd.DataFrame.

fn = "https://raw.githubusercontent.com/PyPSA/powerplantmatching/master/powerplants.csv"
ppl = pd.read_csv(fn, index_col=0)

This dataset includes coordinates (latitude and longitude), which allows us to plot the location and capacity of all power plants in a scatter plot:

ppl.plot.scatter('lon', 'lat', s=ppl.Capacity/1e3)
<AxesSubplot: xlabel='lon', ylabel='lat'>
_images/04-workshop-geopandas_10_1.png

However, this graphs misses some geographic reference point, we’d normally expect for a map like shorelines, country borders etc.

Geopandas - a Pandas extension for geospatial data#

Geopandas extends pandas by adding support for geospatial data.

The core data structure in GeoPandas is the geopandas.GeoDataFrame, a subclass of pandas.DataFrame, that can store geometry columns and perform spatial operations.

Note

Documentation for this package is available at https://geopandas.org/en/stable/.

Typical geometries are points, lines, and polygons. They come from another library called shapely.

First, we need to import the geopandas package. The conventional alias is gpd:

import geopandas as gpd

We can convert the latitude and longitude values given in the dataset to formal geometries (to be exact: shapely.Point objects but we won’t go into detail regarding this) using the gpd.points_from_xy() function, and use this to gpd.GeoDataFrame. We should also specify a so-called coordinate reference system (CRS). The code ‘4326’ means latitude and longitude values.

geometry = gpd.points_from_xy(ppl['lon'], ppl['lat'])
gdf = gpd.GeoDataFrame(ppl, geometry=geometry, crs=4326)

Now, the gdf looks like this:

gdf.head(3)
Name Fueltype Technology Set Country Capacity Efficiency DateIn DateRetrofit DateOut lat lon Duration Volume_Mm3 DamHeight_m StorageCapacity_MWh EIC projectID geometry
id
0 Doel Nuclear Steam Turbine PP Belgium 2911.0 NaN 1975.0 NaN 2022.0 51.32481 4.25889 NaN 0.0 0.0 0.0 {'22WDOELX3000078D', '22WDOELX2000077N', '22WD... {'ENTSOE': {'22WDOELX3000078D', '22WDOELX20000... POINT (4.25889 51.32481)
1 Sarrans Hydro Reservoir Store France 183.0 NaN 1932.0 NaN NaN 44.82942 2.74042 NaN 0.0 0.0 0.0 {'17W100P100P02934'} {'ENTSOE': {'17W100P100P02934'}, 'OPSD': {'OEU... POINT (2.74042 44.82942)
2 Pragneres Hydro Reservoir Store France 189.2 NaN 1953.0 NaN NaN 42.82110 0.01033 NaN 0.0 0.0 0.0 {'17W100P100P02918'} {'ENTSOE': {'17W100P100P02918'}, 'OPSD': {'OEU... POINT (0.01033 42.82110)

With the additional geometry columns, it is now even easier to plot the geographic data:

gdf.plot(
    column='Fueltype',
    markersize=gdf.Capacity/1e2,
)
<AxesSubplot: >
_images/04-workshop-geopandas_24_1.png

We can also start up an interactive map to explore the geodata in more detail:

gdf.explore(column='Fueltype')
Make this Notebook Trusted to load map: File -> Trust Notebook

Map Projections with Cartopy#

Cartopy is a Python package designed for geospatial data processing and has exposed an interface to enable easy map creation using matplotlib.

The Earth is a globe, but we present maps usually on two-dimensional surfaces. Hence, we typically need to project data points onto flat surfaces (e.g. screens, paper). However, we will always loose some information in doing so.

A map projection is:

a systematic transformation of the latitudes and longitudes of locations from the surface of a sphere or an ellipsoid into locations on a plane. Wikipedia: Map projection.

Different projections preserve different metric properties. As a result, converting geodata from one projection to another is a common exercise in geographic data science.

  • conformal projections preserve angles/directions (e.g. Mercator projection)

  • equal-area projections preserve area measure (e.g. Mollweide)

  • equidistant projections preserve distances between points (e.g. Plate carrée)

  • compromise projections seek to strike a balance between distortions (e.g. Robinson)

If you like the “Orange-as-Earth” analogy for projections, checkout this numberphile video by Hannah Fry.

Note

Documentation for this package is available at https://scitools.org.uk/cartopy/docs/latest/.

First, we need to import the relevant parts of the cartopy package:

import cartopy
import cartopy.crs as ccrs

Let’s draw a first map with cartopy outlining the global coastlines in the so-called plate carrée projection (equirectangular projection):

ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7feae46f1850>
Error in callback <function _draw_all_if_interactive at 0x7feb069ceaf0> (for post_execute):
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/pyplot.py:119, in _draw_all_if_interactive()
    117 def _draw_all_if_interactive():
    118     if matplotlib.is_interactive():
--> 119         draw_all()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backend_bases.py:2054, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2052 if not self._is_idle_drawing:
   2053     with self._idle_draw_cntx():
-> 2054         self.draw(*args, **kwargs)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:408, in FigureCanvasAgg.draw(self)
    404 # Acquire a lock on the shared font cache.
    405 with RendererAgg.lock, \
    406      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    407       else nullcontext()):
--> 408     self.figure.draw(self.renderer)
    409     # A GUI class may be need to update a window using this draw, so
    410     # don't forget to call the superclass.
    411     super().draw()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/figure.py:3074, in Figure.draw(self, renderer)
   3071         # ValueError can occur when resizing a window.
   3073 self.patch.draw(renderer)
-> 3074 mimage._draw_list_compositing_images(
   3075     renderer, self, artists, self.suppressComposite)
   3077 for sfig in self.subfigs:
   3078     sfig.draw(renderer)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:543, in GeoAxes.draw(self, renderer, **kwargs)
    535 """
    536 Extend the standard behaviour of :func:`matplotlib.axes.Axes.draw`.
    537 
   (...)
    540 been set.
    541 """
    542 # Shared processing steps
--> 543 self._draw_preprocess(renderer)
    545 # XXX This interface needs a tidy up:
    546 #       image drawing on pan/zoom;
    547 #       caching the resulting image;
    548 #       buffering the result by 10%...;
    549 if not self._done_img_factory:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes._draw_preprocess(self, renderer)
    506 # If data has been added (i.e. autoscale hasn't been turned off)
    507 # then we should autoscale the view.
    508 if self.get_autoscale_on() and self.ignore_existing_data_limits:
--> 509     self.autoscale_view()
    511 # Adjust location of background patch so that new gridlines below are
    512 # clipped correctly.
    513 self.patch._adjust_location()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:946, in GeoAxes.autoscale_view(self, tight, scalex, scaley)
    943 matplotlib.axes.Axes.autoscale_view(self, tight=tight,
    944                                     scalex=scalex, scaley=scaley)
    945 # Limit the resulting bounds to valid area.
--> 946 if scalex and self._autoscaleXon:
    947     bounds = self.get_xbound()
    948     self.set_xbound(max(bounds[0], self.projection.x_limits[0]),
    949                     min(bounds[1], self.projection.x_limits[1]))

AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/IPython/core/formatters.py:339, in BaseFormatter.__call__(self, obj)
    337     pass
    338 else:
--> 339     return printer(obj)
    340 # Finally look for special method names
    341 method = get_real_method(obj, self.print_method)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/IPython/core/pylabtools.py:151, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    148     from matplotlib.backend_bases import FigureCanvasBase
    149     FigureCanvasBase(fig)
--> 151 fig.canvas.print_figure(bytes_io, **kw)
    152 data = bytes_io.getvalue()
    153 if fmt == 'svg':

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backend_bases.py:2314, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2308     renderer = _get_renderer(
   2309         self.figure,
   2310         functools.partial(
   2311             print_method, orientation=orientation)
   2312     )
   2313     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2314         self.figure.draw(renderer)
   2316 if bbox_inches:
   2317     if bbox_inches == "tight":

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/figure.py:3074, in Figure.draw(self, renderer)
   3071         # ValueError can occur when resizing a window.
   3073 self.patch.draw(renderer)
-> 3074 mimage._draw_list_compositing_images(
   3075     renderer, self, artists, self.suppressComposite)
   3077 for sfig in self.subfigs:
   3078     sfig.draw(renderer)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:543, in GeoAxes.draw(self, renderer, **kwargs)
    535 """
    536 Extend the standard behaviour of :func:`matplotlib.axes.Axes.draw`.
    537 
   (...)
    540 been set.
    541 """
    542 # Shared processing steps
--> 543 self._draw_preprocess(renderer)
    545 # XXX This interface needs a tidy up:
    546 #       image drawing on pan/zoom;
    547 #       caching the resulting image;
    548 #       buffering the result by 10%...;
    549 if not self._done_img_factory:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes._draw_preprocess(self, renderer)
    506 # If data has been added (i.e. autoscale hasn't been turned off)
    507 # then we should autoscale the view.
    508 if self.get_autoscale_on() and self.ignore_existing_data_limits:
--> 509     self.autoscale_view()
    511 # Adjust location of background patch so that new gridlines below are
    512 # clipped correctly.
    513 self.patch._adjust_location()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:946, in GeoAxes.autoscale_view(self, tight, scalex, scaley)
    943 matplotlib.axes.Axes.autoscale_view(self, tight=tight,
    944                                     scalex=scalex, scaley=scaley)
    945 # Limit the resulting bounds to valid area.
--> 946 if scalex and self._autoscaleXon:
    947     bounds = self.get_xbound()
    948     self.set_xbound(max(bounds[0], self.projection.x_limits[0]),
    949                     min(bounds[1], self.projection.x_limits[1]))

AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
<Figure size 640x480 with 1 Axes>

A list of the available projections can be found on the Cartopy projection list page.

ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [12], line 2
      1 ax = plt.axes(projection=ccrs.Mollweide())
----> 2 ax.stock_img()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:1064, in GeoAxes.stock_img(self, name)
   1059     source_proj = ccrs.PlateCarree()
   1060     fname = os.path.join(config["repo_data_dir"],
   1061                          'raster', 'natural_earth',
   1062                          '50-natural-earth-1-downsampled.png')
-> 1064     return self.imshow(imread(fname), origin='upper',
   1065                        transform=source_proj,
   1066                        extent=[-180, 180, -90, 90])
   1067 else:
   1068     raise ValueError('Unknown stock image %r.' % name)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:318, in _add_transform.<locals>.wrapper(self, *args, **kwargs)
    313     raise ValueError(f'Invalid transform: Spherical {func.__name__} '
    314                      'is not supported - consider using '
    315                      'PlateCarree/RotatedPole.')
    317 kwargs['transform'] = transform
--> 318 return func(self, *args, **kwargs)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:1368, in GeoAxes.imshow(self, img, *args, **kwargs)
   1364 if not isinstance(transform, ccrs.Projection):
   1365     raise ValueError('Expected a projection subclass. Cannot '
   1366                      'handle a %s in imshow.' % type(transform))
-> 1368 target_extent = self.get_extent(self.projection)
   1369 regrid_shape = kwargs.pop('regrid_shape', 750)
   1370 regrid_shape = self._regrid_shape_aspect(regrid_shape,
   1371                                          target_extent)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:814, in GeoAxes.get_extent(self, crs)
    805 def get_extent(self, crs=None):
    806     """
    807     Get the extent (x0, x1, y0, y1) of the map in the given coordinate
    808     system.
   (...)
    812 
    813     """
--> 814     p = self._get_extent_geom(crs)
    815     r = p.bounds
    816     x1, y1, x2, y2 = r

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:821, in GeoAxes._get_extent_geom(self, crs)
    819 def _get_extent_geom(self, crs=None):
    820     # Perform the calculations for get_extent(), which just repackages it.
--> 821     with self.hold_limits():
    822         if self.get_autoscale_on():
    823             self.autoscale_view()

File ~/micromamba-root/envs/esm/lib/python3.9/contextlib.py:117, in _GeneratorContextManager.__enter__(self)
    115 del self.args, self.kwds, self.func
    116 try:
--> 117     return next(self.gen)
    118 except StopIteration:
    119     raise RuntimeError("generator didn't yield") from None

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:491, in GeoAxes.hold_limits(self, hold)
    488 data_lim = self.dataLim.frozen().get_points()
    489 view_lim = self.viewLim.frozen().get_points()
    490 other = (self.ignore_existing_data_limits,
--> 491          self._autoscaleXon, self._autoscaleYon)
    492 try:
    493     yield

AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
Error in callback <function _draw_all_if_interactive at 0x7feb069ceaf0> (for post_execute):
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/pyplot.py:119, in _draw_all_if_interactive()
    117 def _draw_all_if_interactive():
    118     if matplotlib.is_interactive():
--> 119         draw_all()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backend_bases.py:2054, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2052 if not self._is_idle_drawing:
   2053     with self._idle_draw_cntx():
-> 2054         self.draw(*args, **kwargs)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:408, in FigureCanvasAgg.draw(self)
    404 # Acquire a lock on the shared font cache.
    405 with RendererAgg.lock, \
    406      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    407       else nullcontext()):
--> 408     self.figure.draw(self.renderer)
    409     # A GUI class may be need to update a window using this draw, so
    410     # don't forget to call the superclass.
    411     super().draw()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/figure.py:3074, in Figure.draw(self, renderer)
   3071         # ValueError can occur when resizing a window.
   3073 self.patch.draw(renderer)
-> 3074 mimage._draw_list_compositing_images(
   3075     renderer, self, artists, self.suppressComposite)
   3077 for sfig in self.subfigs:
   3078     sfig.draw(renderer)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:543, in GeoAxes.draw(self, renderer, **kwargs)
    535 """
    536 Extend the standard behaviour of :func:`matplotlib.axes.Axes.draw`.
    537 
   (...)
    540 been set.
    541 """
    542 # Shared processing steps
--> 543 self._draw_preprocess(renderer)
    545 # XXX This interface needs a tidy up:
    546 #       image drawing on pan/zoom;
    547 #       caching the resulting image;
    548 #       buffering the result by 10%...;
    549 if not self._done_img_factory:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes._draw_preprocess(self, renderer)
    506 # If data has been added (i.e. autoscale hasn't been turned off)
    507 # then we should autoscale the view.
    508 if self.get_autoscale_on() and self.ignore_existing_data_limits:
--> 509     self.autoscale_view()
    511 # Adjust location of background patch so that new gridlines below are
    512 # clipped correctly.
    513 self.patch._adjust_location()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:946, in GeoAxes.autoscale_view(self, tight, scalex, scaley)
    943 matplotlib.axes.Axes.autoscale_view(self, tight=tight,
    944                                     scalex=scalex, scaley=scaley)
    945 # Limit the resulting bounds to valid area.
--> 946 if scalex and self._autoscaleXon:
    947     bounds = self.get_xbound()
    948     self.set_xbound(max(bounds[0], self.projection.x_limits[0]),
    949                     min(bounds[1], self.projection.x_limits[1]))

AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/IPython/core/formatters.py:339, in BaseFormatter.__call__(self, obj)
    337     pass
    338 else:
--> 339     return printer(obj)
    340 # Finally look for special method names
    341 method = get_real_method(obj, self.print_method)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/IPython/core/pylabtools.py:151, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    148     from matplotlib.backend_bases import FigureCanvasBase
    149     FigureCanvasBase(fig)
--> 151 fig.canvas.print_figure(bytes_io, **kw)
    152 data = bytes_io.getvalue()
    153 if fmt == 'svg':

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/backend_bases.py:2314, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2308     renderer = _get_renderer(
   2309         self.figure,
   2310         functools.partial(
   2311             print_method, orientation=orientation)
   2312     )
   2313     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2314         self.figure.draw(renderer)
   2316 if bbox_inches:
   2317     if bbox_inches == "tight":

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:74, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     72 @wraps(draw)
     73 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 74     result = draw(artist, renderer, *args, **kwargs)
     75     if renderer._rasterizing:
     76         renderer.stop_rasterizing()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/figure.py:3074, in Figure.draw(self, renderer)
   3071         # ValueError can occur when resizing a window.
   3073 self.patch.draw(renderer)
-> 3074 mimage._draw_list_compositing_images(
   3075     renderer, self, artists, self.suppressComposite)
   3077 for sfig in self.subfigs:
   3078     sfig.draw(renderer)

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/image.py:131, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129 if not_composite or not has_images:
    130     for a in artists:
--> 131         a.draw(renderer)
    132 else:
    133     # Composite any adjacent images together
    134     image_group = []

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/matplotlib/artist.py:51, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     48     if artist.get_agg_filter() is not None:
     49         renderer.start_filter()
---> 51     return draw(artist, renderer)
     52 finally:
     53     if artist.get_agg_filter() is not None:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:543, in GeoAxes.draw(self, renderer, **kwargs)
    535 """
    536 Extend the standard behaviour of :func:`matplotlib.axes.Axes.draw`.
    537 
   (...)
    540 been set.
    541 """
    542 # Shared processing steps
--> 543 self._draw_preprocess(renderer)
    545 # XXX This interface needs a tidy up:
    546 #       image drawing on pan/zoom;
    547 #       caching the resulting image;
    548 #       buffering the result by 10%...;
    549 if not self._done_img_factory:

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:509, in GeoAxes._draw_preprocess(self, renderer)
    506 # If data has been added (i.e. autoscale hasn't been turned off)
    507 # then we should autoscale the view.
    508 if self.get_autoscale_on() and self.ignore_existing_data_limits:
--> 509     self.autoscale_view()
    511 # Adjust location of background patch so that new gridlines below are
    512 # clipped correctly.
    513 self.patch._adjust_location()

File ~/micromamba-root/envs/esm/lib/python3.9/site-packages/cartopy/mpl/geoaxes.py:946, in GeoAxes.autoscale_view(self, tight, scalex, scaley)
    943 matplotlib.axes.Axes.autoscale_view(self, tight=tight,
    944                                     scalex=scalex, scaley=scaley)
    945 # Limit the resulting bounds to valid area.
--> 946 if scalex and self._autoscaleXon:
    947     bounds = self.get_xbound()
    948     self.set_xbound(max(bounds[0], self.projection.x_limits[0]),
    949                     min(bounds[1], self.projection.x_limits[1]))

AttributeError: 'GeoAxesSubplot' object has no attribute '_autoscaleXon'
<Figure size 640x480 with 1 Axes>

We can combine the functionality of cartopy with geopandas plots:

fig = plt.figure(figsize=(7,7))

ax = plt.axes(projection=ccrs.PlateCarree())

gdf.plot(
    ax=ax,
    column='Fueltype',
    markersize=gdf.Capacity/1e2,
)
<GeoAxesSubplot:>
_images/04-workshop-geopandas_39_1.png

We can add further geographic features to this map for better orientation.

For instance, we can add the coastlines…

ax.coastlines()
fig
_images/04-workshop-geopandas_41_0.png

… country borders …

ax.add_feature(cartopy.feature.BORDERS, color='grey', linewidth=0.5)
fig
_images/04-workshop-geopandas_43_0.png

… colour in the ocean in blue …

ax.add_feature(cartopy.feature.OCEAN, color='azure')
fig
_images/04-workshop-geopandas_45_0.png

…and color in the land area in yellow …

ax.add_feature(cartopy.feature.LAND, color='cornsilk')
fig
_images/04-workshop-geopandas_47_0.png

Geopandas will automatically calculate sensible bounds for the plot given the geographic data. But we can also manually zoom in or out by setting the spatial extent with the .set_extent() method:

ax.set_extent([5, 16, 47, 55])
fig
_images/04-workshop-geopandas_49_0.png

Reprojecting a GeoDataFrame#

In geopandas, we can use the function .to_crs() to convert a GeoDataFrame to a desired coordinate reference system. In this particular case, we use the proj4_init string of an initialised cartopy projection to reproject our power plant GeoDataFrame.

fig = plt.figure(figsize=(7,7))

crs = ccrs.AlbersEqualArea()

ax = plt.axes(projection=crs)

gdf.to_crs(crs.proj4_init).plot(
    ax=ax,
    column='Fueltype',
    markersize=gdf.Capacity/1e2,
)

ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fbd9d28a070>
_images/04-workshop-geopandas_52_1.png

Reading and Writing Files with geopandas#

In the following example, we’ll load a dataset containing the NUTS regions:

Nomenclature of Territorial Units for Statistics or NUTS (French: Nomenclature des unités territoriales statistiques) is a geocode standard for referencing the subdivisions of countries for statistical purposes.

Our ultimate goal for this part of the tutorial is to map the power plant capacities to the NUTS-1 region they belong to.

Common filetypes for vector-based geospatial datasets are GeoPackage (.gpkg), GeoJSON (.geojson), File Geodatabase (.gdb), or Shapefiles (.shp).

In geopandas we can use the gpd.read_file() function to read such files. So let’s start:

nuts = gpd.read_file("../../data/nuts/NUTS_RG_10M_2021_4326.geojson")
nuts.head(3)
id NUTS_ID LEVL_CODE CNTR_CODE NAME_LATN NUTS_NAME MOUNT_TYPE URBN_TYPE COAST_TYPE FID geometry
0 BG423 BG423 3 BG Pazardzhik Пазарджик 3.0 2 3 BG423 POLYGON ((24.42101 42.55306, 24.41032 42.46950...
1 BG424 BG424 3 BG Smolyan Смолян 3.0 3 3 BG424 POLYGON ((25.07422 41.79348, 25.05851 41.75177...
2 BG425 BG425 3 BG Kardzhali Кърджали 3.0 3 3 BG425 POLYGON ((25.94863 41.32034, 25.90644 41.30757...

It is good practice to set an index. You can use .set_index() for that:

nuts = nuts.set_index('id')

We can also check out the geometries in the dataset with .geometry:

nuts.geometry
id
BG423    POLYGON ((24.42101 42.55306, 24.41032 42.46950...
BG424    POLYGON ((25.07422 41.79348, 25.05851 41.75177...
BG425    POLYGON ((25.94863 41.32034, 25.90644 41.30757...
CH011    MULTIPOLYGON (((6.86623 46.90929, 6.89621 46.9...
CH012    POLYGON ((8.47767 46.52760, 8.39953 46.48872, ...
                               ...                        
LV       POLYGON ((27.35158 57.51824, 27.54521 57.53444...
ME       POLYGON ((20.06394 43.00682, 20.32958 42.91149...
MK       POLYGON ((22.36021 42.31116, 22.51041 42.15516...
SK0      POLYGON ((19.88393 49.20418, 19.96275 49.23031...
IT       MULTIPOLYGON (((12.47792 46.67984, 12.69064 46...
Name: geometry, Length: 2010, dtype: geometry

With .crs we can check in which coordinate reference system the data is given:

nuts.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
nuts.total_bounds
array([-63.08825, -21.38917,  55.83616,  80.76427])

Let’s filter by NUTS-1 level…

nuts1 = nuts.query("LEVL_CODE == 1")

… and explore what kind of geometries we have in the dataset …

nuts1.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook

To write a GeoDataFrame back to file use GeoDataFrame.to_file(). The file format is inferred from the file ending.

nuts1.to_file("NUTS1.geojson")

Calculating the areas and buffers#

The first thing we need to do to calculate area or buffers is to reproject the GeoDataFrame to an equal-area projection (here: EPSG:3035):

nuts1 = nuts1.to_crs(3035)

The area can be accessed via .area and is given in m² (after projection). Let’s convert to km²:

area = nuts1.area / 1e6
area
id
AT1    23545.286205
AT2    25894.953057
EL4    17388.679384
EE0    45315.713593
EL3     3799.676547
           ...     
PL7    29846.398582
PL8    63217.536546
PL9    35563.812826
RO2    72545.290405
SK0    49008.115415
Length: 125, dtype: float64
nuts1.explore(column=area, vmax=1e5)
Make this Notebook Trusted to load map: File -> Trust Notebook

We can also build a buffer of 1km around each geometry using .buffer():

nuts1.buffer(1000).explore()
Make this Notebook Trusted to load map: File -> Trust Notebook

Joining spatial datasets#

Multiple GeoDataFrames can be combined via spatial joins.

Observations from two datasets are combined with the .sjoin() function based on their spatial relationship to one another (e.g. whether they are intersecting or overlapping). You can read more about the specific options here.

Let’s first reproject the gdf object to the same CRS as nuts1:

gdf = gdf.to_crs(3035)

Then, let’s have a look at both datasets at once. We want to find out which points (representing power plants) lie within which shape (representing NUTS regions).

fig = plt.figure(figsize=(7,7))

ax = plt.axes(projection=ccrs.epsg(3035))

nuts1.plot(
    ax=ax,
    edgecolor='black',
    facecolor='lightgrey'
)

gdf.to_crs(3035).plot(
    ax=ax,
    column='Fueltype',
    markersize=gdf.Capacity/20,
    legend=True
)

ax.set_extent([5, 19, 47, 55])
_images/04-workshop-geopandas_85_0.png

We can now apply the .sjoin function to look for which power plants lie within which NUTS1 region. By default, sjoin looks for intersections and keeps the geometries of the left GeoDataFrame.

joined = gdf.sjoin(nuts1)

If we look at this new GeoDataFrame, we now have additional columns from the NUTS1 data:

joined.head(3)
Name Fueltype Technology Set Country Capacity Efficiency DateIn DateRetrofit DateOut ... index_right NUTS_ID LEVL_CODE CNTR_CODE NAME_LATN NUTS_NAME MOUNT_TYPE URBN_TYPE COAST_TYPE FID
id
0 Doel Nuclear Steam Turbine PP Belgium 2911.0 NaN 1975.0 NaN 2022.0 ... BE2 BE2 1 BE Vlaams Gewest Vlaams Gewest 0.0 0 0 BE2
170 Drogenbos Tgv Natural Gas CCGT PP Belgium 465.0 NaN NaN NaN NaN ... BE2 BE2 1 BE Vlaams Gewest Vlaams Gewest 0.0 0 0 BE2
172 Rodenhuize Bioenergy Steam Turbine PP Belgium 268.0 NaN NaN NaN NaN ... BE2 BE2 1 BE Vlaams Gewest Vlaams Gewest 0.0 0 0 BE2

3 rows × 29 columns

We can now use these new columns to group the capacities (and convert to a suitable unit):

cap = joined.groupby("NUTS_ID").Capacity.sum() / 1000 # GW

Let’s quickly check if all NUTS1 regions have power plants:

nuts1.index.difference(cap.index)
Index(['CY0', 'ES7', 'FI2', 'FRY', 'IS0', 'LI0', 'MK0', 'MT0', 'PT2', 'PT3',
       'TR1', 'TR2', 'TR3', 'TR4', 'TR5', 'TR6', 'TR7', 'TR8', 'TR9', 'TRA',
       'TRB', 'TRC'],
      dtype='object')

This is not the case. Then it is good practice to reindex the series to include all NUTS1 regions, even if this leads to some NaN values.

cap = cap.reindex(nuts1.index)
cap
id
AT1    4.382100
AT2    4.089500
EL4    0.427967
EE0    1.369000
EL3    0.022142
         ...   
PL7    6.379902
PL8    0.685651
PL9    5.230788
RO2    2.036566
SK0    5.567220
Name: Capacity, Length: 125, dtype: float64

Finally, we can plot the total generation capacity per NUTS1 region on a map.

nuts1.plot(figsize=(7,7), column=cap, legend=True)
<AxesSubplot:>
_images/04-workshop-geopandas_98_1.png

This concludes the geopandas and cartopy tutorial.

Exercises#